### Staggered Diff-in-diff

## Reinstall version 0.8.1 for reproducibility
#library(versions)
#vers <- "0.8.1"
#install.versions('fixest', vers)

rm(list = ls())
setwd("C:/Users/natos/Box/Trejo & Skigin/Journalists/data")
library(tidyverse)
library(fixest)
library(haven)
here::here()

d <- as_tibble(read.dta13("data-state-level.dta"))

# year treated, control group set to 1000
data.state$year_treat = ifelse(data.state$statecode == 2, 2007,  # BajaCalif
                               ifelse(data.state$statecode == 8, 2008, #Chihuahua
                                      # Durango
                                      ifelse(data.state$statecode == 10, 2008,
                                             #Guerrero
                                             ifelse(data.state$statecode == 12, 2007,
                                                    # Michoacan
                                                    ifelse(data.state$statecode == 16, 2007,
                                                           #Nuevo Leon
                                                           ifelse(data.state$statecode == 19, 2008,
                                                                  # Sinaloa
                                                                  ifelse(data.state$statecode == 25, 2008,
                                                                         # Tamaulipas
                                                                         ifelse(data.state$statecode == 28, 2008,
                                                                                ifelse(data.state$statecode == 30, 2011, 1000)))))))))



# time to treatment
data.state$time_to_treat = data.state$year-data.state$year_treat
data.state$time_to_treat = ifelse(data.state$year_treat ==1000,-1000,data.state$time_to_treat)


# with cohort x time to treatment dummies
res_cohort = feols(count ~ i(time_to_treat, f2 = statecode, drop = c(-1000)) | 
                     statecode + year, data.state)



# SA method: we aggregate the effects for each period
agg_coef = aggregate(res_cohort, "(ti.*at)::(-?[[:digit:]])")[-1,]
x = c(-9:-2, 0:7)

pe.minus1 = aggregate(res_cohort, "(ti.*at)::(-?[[:digit:]])")[1,1]
ci_lo.minus1 = aggregate(res_cohort, "(ti.*at)::(-?[[:digit:]])")[1,1] -
  1.96 * aggregate(res_cohort, "(ti.*at)::(-?[[:digit:]])")[1,2]
ci_up.minus1 = aggregate(res_cohort, "(ti.*at)::(-?[[:digit:]])")[1,1] +
  1.96 * aggregate(res_cohort, "(ti.*at)::(-?[[:digit:]])")[1,2]

## Figure 8
pdf("Figure8.pdf",
    width = 17,
    height = 8.5)
png("../figures/submission/Figure8.png",
    width = 465, height = 225, units='mm', res = 600)
par(mar = c(9, 6, 1.5, .5))
plot(x, agg_coef[, 1], pch = 17, col = 4, ylim = c(-1,2),
     xlab = "", ylab = "", cex = 2, xaxt="n", cex.axis=1.5)
axis(1, at = seq(-9, 7, by = 1), las=1,cex.axis=1.5)
abline(h=seq(-1,2,.5), lty=2, col="gray")
abline(h=0, lty=2)
abline(v=-9:8, lty=2, col="gray")
abline(v=-.2, lty=2, lwd=2)
ci_low = agg_coef[, 1] - 1.96 * agg_coef[, 2]
ci_up = agg_coef[, 1] + 1.96 * agg_coef[, 2]
segments(x0 = x, y0 = ci_low, x1 = x, y1 = ci_up, col = 4, cex = 2, lwd=3)
# Manually add beta and CI for x = -1
points(-1, pe.minus1, pch = 17, col = 4, cex = 2)
segments(x0 = -1, y0 = ci_lo.minus1, x1 = -1, y1 = ci_up.minus1, col = 4,
         cex = 2, lwd=3)
title(ylab = expression(italic("Effect of Militarization")), line = 2.7,
      cex.lab = 2)
title(xlab = expression(italic("Time to treatment")), line = 3.5, cex.lab = 2)
text(x = .85, y = -.3, "First year of \nmilitarization", cex = 1.8)
dev.off()
